Bolker et al. ( ???): Generalized linear mixed models: a practical guide for ecology and evolution

Very Important:

• Understand the nature of the problems we are trying to solve
• Understand the value in solving them

Moderately Important:

• Understand some of the technical aspects of multilevel model fitting and interpretation

Minimally Important:

• Understand the math involved

1. Normality (check pooled residuals)
2. Homogeneity (check pooled residuals, transformation?)
3. Fixed explanatory variables (plan sampling)
4. Independence (observations often correlated)
5. Correct model

Our tests (so far) assume independence of samples.

Lots of data have structure at additional levels beyond the main variable(s) of interest

• Error in these levels is correlated (making the data non-independent)

We need a method to model that correlation structure.

Explicitly by design:

• Grouping (leaves within plants within treatments within fields)
• Repeated measurements of the same individual (time or space)

Implicitly as a side effect:

• Non-equal sample sizes, missing data
• Non-equal within-group variance
• Phylogenetic structure ( ???)
• Population structure and genetic relatedness ( ???)

Multilevel / hierarchical / mixed / repeated measures models

• All refer to the same class of models, which contain both fixed and random variables.
• Some will further subdivide and draw fine distinctions between models.

Fixed effects:

• Explanatory variables of primary interest, which would be used again if the experiment were repeated
• Sex, treatment groups, etc.

Random effects

• Variables that are a "random" selection from a larger set of variables that allow for heterogeneity, correlation, phylogenetic or genetic structure, real random variation (\(\epsilon\)).
• Subject, block, etc.
• "Random" does not necessarily mean "not of interest."

Also see: http://andrewgelman.com/2005/01/25/why_i_dont_use/

Units are grouped at different levels.

• Grouping takes different forms

\[Y = \beta_0 + \beta_1 X_1 + ... + \beta_n X_n + \mbox{Random}\]

Where \(\mbox{Random}\) can be:

• Random noise (\(\epsilon\))
• Heterogeneity of variance
• Nesting structure
• Temporal correlation (autocorrelation)
• Spatial correlation
• Phylogenetic correlation
• Genetic correlation

• Maneuver (hover vs. clockwise, hover vs. counterclockwise)
  • Bird (4)
    • Trial (5)
      • Wingbeat (mean of 16 \(\rightarrow\) 1)

1. Maneuver is the variable of interest ("fixed")
2. Birds are a representative sample of all birds ("random")
3. Trials are repeated measures of each bird ("random")
4. Wingbeats are the observation level

1. Calculate means of 16 wingbeats for each bird for each trial
   • Loss of information, but…
2. Five trials per bird
3. Multilevel ANOVA
   • Bird and trial nested within bird as random effects
   • Maneuver as the fixed effect

Each bird gets its own mean value for a maneuver (across 5 trials). Maneuver is analyzed as the mean difference between each maneuver.

Conceptual complexity

• You must know your model
• Difficult to visualize
• Do my results make sense?

Computational complexity

• Many methods fit by likelihood
  • Singularities abound, differentiation doesn't always work
• MCMC can be slow, difficult to tune
• Bayesian requires defining priors

Among the options for R packages:

• nlme: Oldest, restricted to a subset of models (linear and non-linear Gaussian), for us the place to start
• lme4: Newer (but same authors as nlme), can handle more complex models with crossed random effects, pedigrees, non-Gaussian response, etc., more flexible distributions of response variables, but harder to test hypotheses
• MCMCglmm, rstan, rjags: Fit Bayesian multilevel models via a variety of MCMC samplers

Intraclass correlation coefficient:

• Repeated measurements for a single specimen
• What is the relative among measurement vs. between specimen variation

library(nlme)
fm_lme <- lme(Fem_Len ~ 1,
              random = ~ 1 | MouseID,
              data = Lengths,
              method = "ML")
head(coef(fm_lme))

##       (Intercept)
## 46001    16.21499
## 46002    17.08474
## 46003    16.77983
## 46004    16.64986
## 46012    16.85481
## 46013    15.99005

Lengths %>% group_by(MouseID) %>% summarize(mean(Fem_Len)) %>% head()

## # A tibble: 6 × 2
##   MouseID `mean(Fem_Len)`
##    <fctr>           <dbl>
## 1   46001          16.215
## 2   46002          17.085
## 3   46003          16.780
## 4   46004          16.650
## 5   46012          16.855
## 6   46013          15.990

(varcomps <- VarCorr(fm_lme))

## MouseID = pdLogChol(1) 
##             Variance     StdDev    
## (Intercept) 0.3151310080 0.56136531
## Residual    0.0001810924 0.01345706

var.among <- as.numeric(varcomps[1, 1])
var.within <- as.numeric(varcomps[2, 1])
var.among / (var.among + var.within)

## [1] 0.9994257

Note: very close but not numerically equivalent to the anova() version. Models are slightly different.

Seasonal patterns in testis mass in the squid Loligo. Testis mass predicted by dorsal mantle length, month, and DML X month.

S <- read_excel("../data/Squid.xlsx")
S <- S %>% mutate(Month = factor(Month))
str(S)

## Classes 'tbl_df', 'tbl' and 'data.frame':    768 obs. of  5 variables:
##  $ Specimen   : num  1017 1034 1070 1070 1019 ...
##  $ Year       : num  1991 1990 1990 1990 1990 ...
##  $ Month      : Factor w/ 12 levels "1","2","3","4",..: 2 9 12 11 8 10 5 7 7 7 ...
##  $ DML        : num  136 144 108 130 121 117 133 105 109 97 ...
##  $ Testis_Mass: num  0.006 0.008 0.008 0.011 0.012 0.012 0.013 0.015 0.017 0.017 ...

# Check for missing data
sum(!complete.cases(S))

## [1] 0

Testis mass modeled by DML, Month, and the DML X Month interaction:

fm1 <- lm(Testis_Mass ~ DML + Month + DML:Month, data = S)

Extract the residuals and fitted values:

S$Res1 <- residuals(fm1)
S$Fitted <- fitted.values(fm1)

What we need:

1. Allow for change in variance with values of DML (here, more spread with larger values)
2. Allow for difference variances per month

Solution:

• Generalized least squares \(\approx\) "Weighted regression"
• nlme::gls()
• Also one solution for phylogenetic least squares (coming soon)

library(nlme)
variance_struct <- varPower(form = ~ DML | Month)

• This structure allows the residual variance to change with DML on a per-Month basis.
• Here a power relationship
  • Other variance structures available

fm2 <- gls(Testis_Mass ~ DML * Month, data = S,
           weights = variance_struct)

plot(fm2)

S$Res2 <- residuals(fm2, type = "normalized")

Data with hierarchical structure:

• i samples within n blocks

We could just take a means of each block and use those data for ANOVA.

• Lose all the information about variability within blocks
• Might increase apparent variance

Collected by National Institute for Coastal and Marine Management (RIKZ)

RIKZ <- read_excel("../data/RIKZ.xlsx")
glimpse(RIKZ)

## Observations: 45
## Variables: 5
## $ Sample   <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16...
## $ Richness <dbl> 11, 10, 13, 11, 10, 8, 9, 8, 19, 17, 6, 1, 4, 3, 3, 1...
## $ Exposure <dbl> 10, 10, 10, 10, 10, 8, 8, 8, 8, 8, 11, 11, 11, 11, 11...
## $ NAP      <dbl> 0.045, -1.036, -1.336, 0.616, -0.684, 1.190, 0.820, 0...
## $ Beach    <dbl> 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4,...

• Richness: Number of species
• NAP: Height above mean tidal level
• Exposure: Index composed of wave action, length of the surf zone, slope, grain size, and the depth of the anaerobic layer

Convert Beach to factor.

RIKZ <- RIKZ %>% mutate(Beach = factor(Beach))

Sites are the level of the observation.

RIKZ %>% group_by(Beach) %>% tally()

## # A tibble: 9 × 2
##    Beach     n
##   <fctr> <int>
## 1      1     5
## 2      2     5
## 3      3     5
## 4      4     5
## 5      5     5
## 6      6     5
## 7      7     5
## 8      8     5
## 9      9     5

5 observations at each of 9 beaches

Each Beach only has 1 value of Exposure.

fm1 <- lm(Richness ~ NAP + Beach, data = RIKZ)
Anova(fm1, type = "III")

## Anova Table (Type III tests)
## 
## Response: Richness
##             Sum Sq Df F value    Pr(>F)    
## (Intercept) 466.37  1 49.8039 3.224e-08 ***
## NAP         230.66  1 24.6322 1.793e-05 ***
## Beach       416.37  8  5.5581 0.0001441 ***
## Residuals   327.74 35                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Maybe we are not interested in knowing the exact nature of the beach effect.

Allow a different mean richness per beach:

fm2 <- lme(Richness ~ NAP, random = ~ 1 | Beach, data = RIKZ)

• Similar to the model we used for intraclass correlation coefficient.

plot(fm2)

Generate predicted values for overall population level estimates and within beach estimates.

# Overall
RIKZ$Pred0 <- fitted(fm2, level = 0)

# Within beach
RIKZ$Pred1 <- fitted(fm2, level = 1)

Allow each Beach to have its own NAP to Richness linear relationship:

fm3 <- lme(Richness ~ NAP, random = ~ NAP | Beach, data = RIKZ)

Now the NAP to Richness relationship can vary for each Beach.

Model formula for random effects:

\[\sim \mbox{Slope}~|~\mbox{Intercept}\]

Generate predicted values for overall population level estimates and within beach estimates.

# Overall
RIKZ$Pred0 <- fitted(fm3, level = 0)

# Within beach
RIKZ$Pred1 <- fitted(fm3, level = 1)

##   (Intercept)       NAP
## 1    8.421057 -3.656268
## 2   12.363497 -5.536814
## 3    3.806643 -1.505716
## 4    3.562422 -1.385960
## 5   11.200150 -5.137323
## 6    4.426279 -1.775675
## 7    4.082944 -1.644329
## 8    5.099895 -2.106852
## 9    6.335471 -2.721314

Sequential measurements of the same individual / plot / organism, etc. will result in within-subject correlation.

• Growth curves
• Long-term ecological monitoring

Distance from the sella turcica to the pterygomaxillary fissure in growing children.

O <- read_excel("../data/Ortho.xlsx")
O <- O %>% mutate(Subject = factor(Subject),
                  Sex = factor(Sex))
glimpse(O)

## Observations: 108
## Variables: 4
## $ Distance <dbl> 26.0, 25.0, 29.0, 31.0, 21.5, 22.5, 23.0, 26.5, 23.0,...
## $ Age      <dbl> 8, 10, 12, 14, 8, 10, 12, 14, 8, 10, 12, 14, 8, 10, 1...
## $ Subject  <fctr> M01, M01, M01, M01, M02, M02, M02, M02, M03, M03, M0...
## $ Sex      <fctr> Male, Male, Male, Male, Male, Male, Male, Male, Male...

O %>% group_by(Sex) %>% tally()

## # A tibble: 2 × 2
##      Sex     n
##   <fctr> <int>
## 1 Female    44
## 2   Male    64

range(O$Age)

## [1]  8 14

Potthoff, R. F. and Roy, S. N. (1964) A generalized multivariate analysis of variance model useful especially for growth curve problems. Biometrika 51: 313-326.

• Main effects of Age and Sex (like an ANCOVA model).
• Random effects allow the Age by Distance relationship (slope) to vary for each Subject, which is also allowed to have its own mean.

fm2 <- lme(Distance ~ Sex + Age, random = ~ Age | Subject,
           method = "ML",
           data = O)

# Overall
O$Pred0 <- fitted(fm2, level = 0)

# Age within Subject
O$Pred1 <- fitted(fm2, level = 1)

Anova(fm2, type = "III")

## Analysis of Deviance Table (Type III tests)
## 
## Response: Distance
##                Chisq Df Pr(>Chisq)    
## (Intercept) 283.5863  1  < 2.2e-16 ***
## Sex           8.6649  1   0.003244 ** 
## Age          89.1482  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fm1 <- lme(Distance ~ Sex + Age, random = ~ 1 | Subject, 
           method = "ML",
           data = O)

# Overall
O$Pred0 <- fitted(fm1, level = 0)

# Age within Subject
O$Pred1 <- fitted(fm1, level = 1)

anova(fm1, fm2)

##     Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## fm1     1  5 444.8565 458.2671 -217.4282                        
## fm2     2  7 446.8352 465.6101 -216.4176 1 vs 2 2.021324   0.364

No.

Watch Lecture 08-4